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Abstract 



Understanding the time evolution of fragmented animal populations and their habitats, connected by migration, is a problem of 
both theoretical and practical interest. This paper presents a method for calculating the time evolution of the habitats' population 
size distribution from a general stochastic dynamic within each habitat, using a deterministic approximation which becomes exact 
for an infinite number of habitats. Fragmented populations are usually thought to be characterized by a separation of time scale 
between, on the one hand, colonization and extinction of habitats and, on the other hand, the local population dynamics within 
each habitat. The analysis in this paper suggests an alternative view: the effective population dynamic stems from a law of large 
numbers, where stochastic fluctuations in population size of single habitats are buffered through the dispersal pool so that the global 
population dynamic remains approximately smooth. For illustration, the deterministic approximation is compared to simulations of 
a stochastic model with density dependent local recruitment and mortality. The article is concluded with a discussion of the general 
implications of the results, and possible extensions of the method. 
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1. Introduction 

I The habitats of many animal populations are geographically 
I ^ . 'divided into numerous small patches, either because of human 
activities or because the animals natural habitats are patchy. 
Understanding the dynamics of such populations is a problem 
of great theoretical and practical interest. Isolated small habi- 
tats are often extinction prone, because of, e.g., inbreeding in 
.co mbination with demographic and environmental stochastic- 
ity JHanskiL \l99^ . When the habitats are connected into a 
network by migration, local populations may still be extinc- 
tion prone, but the whole population can be sustained because 
empty habitats are recolonized by migrants from surrounding 
occupied habitats. When the migration rate is high, also the lo- 
J> cal populations are stabilized by an inflow of immigrants (the 
rescue effect). 

rS There are two main approaches to modeling such popula- 
tions: First, simplified models that may be treated analytically, 
and, second, more detailed models where simulations are nec- 
essary. In the framework of LevinI ( 1974 ). the full population 
dynamic is simplified by treating each patch as either occupied 
or empty. The rates of patch colonization and extinction are 
expressed as functions of the total fraction of occupied patches 
in the population, or as functions of the number of neighbour- 
ing patches in spatially explicit models. This reduction is gen- 
erally motivated by a separation of time scale between colo- 
nization and extincti on of habitats and the population dynamics 
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within the habitats dUevinl 1 19741: iHanski & GvUenberj, Il993 
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Lande et al. L 119981) . Many of the mechanisms that can be ob- 
served in natural populations (e.g. AUee and rescue effects) 
have correspondences in Levin's framework, and the qualita- 
tive behaviour of such model s is well understood (se e , e.g. , 
i Hanskil 1 1998^ lEtiennel [2OO0I: iHarding & McNamaral 12002 : 
'zhou & Wang 2004*: 'Zhouefa/.L 12004 iTaylor & Hastings . 
J005; Martcheya & Bolker, 2007) . More detailed models in- 
corporate e.g. the effect of correlations between neighbour- 
ing patc hes, leadi n g to non-uniform patter ns of colonized 



patches ( Hui &Lil 12004 iRov et al. I 120081) , and of demo 



graphic stochasticity within the local habitats ( Keelingl 2002). 
The study of such models generally requires computer simula- 
tions (but see Hui & Li, 2004, for an exception). 

Models that describe the metapopulation dynamic as a func- 
tion of the fraction of occupied patches have no explicit con- 
nection to the population dynamic within each patch, i.e. 
it is not apparent how individual births, deaths, or migra- 
tion events translate into effective colonization a nd extinction 
rates. Starting with iHanski & GyllenbergI jl993h . several au- 
thors have attempted to bridge the gap between detailed lo- 
cal population dynamic and the dynamic at the overall popu- 



lation level (e.g.'Prechsler & Wissel'.*^1 997l : lLande et al. 
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Nachman, 2000; Lopez & Pfister, 200l!). This paper contributes 
to these efforts by deriving, from first principles, how the global 
stochastic dynamic leads to an effective deterministic time evo- 
lution of the habitats' population size distribution. This gives 
additional insight in the dynamics of fragmented populations. 
The only critical assumption in this derivation is that the num- 
ber of habitats is large enough. Especially, the method can be 
used for any local stochastic population dynamic. 

The rest of the paper is organized as follows: In section |2] 
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the stochastic population model is presented; in section |3] the 
method is developed from the full stochastic dynamic of the 
population; in section |4] the method is compared to stochastic 
simulations of a model with local density dependent recruit- 
ment and mortality. Finally, section|5]discusses the general im- 
plications of the results, and possible extensions of the method. 

2. Population model 

Consider a population with habitats (patches), all equal. 
In each patch, the population changes according to a Poisson 
process, where the rate of transition from / individuals to /' in- 
dividuals is given by Tji. As an important special case, we have 
the standard birth-death dynamic, where individuals are born 
and die according to Poisson processes with rates B, and D,, re- 
spectively, where / is the number of individuals in the patch. In 
the general model, this corresponds to taking 
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Figure 1 : Illustration of the population dynamic of a patchy population. In each 
habitat, there are 4- individuals in habitat k, and M emigrants in the dispersal 
pool. 



Bj when /' = / - 1 
Dj when j = i + 1, 
else. 



(1) 



In addition to the local population dynamic, the number of 
individuals in each habitat can change because some indi- 
viduals emigrate from the habitat to other patches, or be- 
cause immigrants arrive f rom other patches. Following 



Hanski & GvllenbergI (Il993h . the migration process is modeled 
as follows: Individuals emigrate from a patch at rate where / 
is the patch population size. If the individuals migrate indepen- 
dently and with constant rate, £, is proportional to /, but more 
complex migration patterns can be incorporated. If for exam- 
ple individuals move to avoid overcrowding, the migration rate 
becomes density dependent. The emigrants enter a common 
dispersal pool, containing the migrants from all patches that 
have not yet reached their target patch. Each migrant stays an 
exponentially distributed time in the pool, with expected value 
1/a, before reaching the target habitat, which is chosen with 
equal probability among all habitats. This process is illustrated 
in Figure [U 

Migration may fail if the individuals die before reaching the 
new habitat; this is modeled by the rate v of dying per unit of 
time during dispersal. In practice, the probability of success- 
ful migration depends on the background mortality of the in- 
dividuals, on the time it spends in the dispersal pool, and on 
additional perils that the individuals may be exposed to during 
dispersal (e.g. increased risk of predation due to lack of cover, 
etc.). Thus, if there are M migrants, (a + v)M individuals leave 
the dispersal pool per unit of time, and the rate of immigration 
to a given patch is / = aM/N. 

3. Approximative dynamic 

Since the local population dynamic is assumed to be the same 
within all habitats, it is sufficient to count the number of habitats 
with a given number of individuals rather than keeping track of 
the number of individuals in each habitat. Let n, denote the 



number of habitats with / inhabitants at a given time. Since 
the total number of patches is constant, we have ~ ^ 

at all times. When the number of individuals change from i to 
j in some habitat, n, decreases with one and nj increases with 
one. From the model in section |2] we can write down a set 
of Master equations for the probability p of observing counts 
no,ni, . . . and M individuals in the dispersal pool (see the ap- 
pendix). In principle, these equations can be integrated ana- 
lytically. In practice, however, because the number of possible 
configurations grows exponentially with the number of patches, 
it is not possible to use the full Master equation except when the 
number of habitats and the maximum population size are both 
small. 

When the number of habitats is large enough that the total 
population size is large, even though the number of individu- 
als per habitat may be small, we can use a different approach. 
In this case, we can expect that the probability p of observing a 
given distribution of habitat population sizes and number of mi- 
grants changes only a little when a single habitat (or a habitat 
and the migration pool) changes its population size. It is im- 
portant here to emphasize that we make no assumptions on the 
size of the changes to the population size of any single habitat, 
only on that there are sufficiently many habitats that they form 
a statistical ensemble. Hence, on the level of the population, 
this leads only to a single increment and decrement in the count 
of the number of patches with a given number of individuals. 
For instance, it is perfectly possible to have big jumps in the 
population number of individual patches in the model without 
violating the smoothness of the distribution of habitat popula- 
tion size. 

It is now convenient to express the frequency of habitat pop- 
ulation sizes in terms of the scaled frequency / - rij/N and 
the patch immigration rate / = aM/N. The probability p of 
observing /o, f\, ■ ■ ■ and / is related to the probability p by 



p(fo,fu...J)^p(Nfo,Nfu 



,NI/a). 



(2) 



In contrast to n, and M, the distribution of fi and / can be ex- 
pected to become approximately independent of for a large 
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number of habitats. A change of a single count n, leads only to 
the change 1 IN in fj. Similarly, a single increment of M leads 
to the change a IN in /. Since for large the probability p is an 
approximately continuous function of /j and /, we can expand 
the full Master equation of p in terms of A^"' to obtain a par- 
tial differential equation for p (see Appendix). Taking only the 
leading term (ignoring terms of order A^"' or higher), we obtain 
a deterministic dynamic for fi and /, given by 
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d,I^-{a + v)I + aY^Ekfk. 



(3) 
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In the limit of infinitely many patches, Eq. (O gives an exact 
description of the dynamic. 

When the migrants spend a short time in the dispersal state 
compared to the normal birth-death dynamic, i.e. when a is 
large, / rapidly adjusts so that dj becomes approximately zero, 
i.e. 
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(4) 



at each moment in time. Inserting this into Eq. ^ yields a 
nonlinear Master equation for the patch dynamic. 
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In this equation, a I {a + v) is the probability that a migrant suc- 
cessfully reaches a new habitat. Depending on the biological 
details of the dispersal, which are not explicitly modeled, v may 
depend on a. The stationary solution of Eq. |5]is similar to the 
result by ,Nach man (2000). He does not, however, derive the 
equations from the full stochastic dynamic. In the next sec- 
tion, I investigate the accuracy of this approximation by com- 
paring the solution of the non-linear Master equation to explicit 
stochastic simulations of a population where the local dynamic 
is governed by simple growth and density dependent mortality. 

4. Simple density dependent dynamics 

As a simple but illustrative example, we consider a popula- 
tion with constant birth rate and migration rate per capita, and 
density dependent mortality. The birth, death, and migration 
rates of a patch with ; individuals are 



Bi = ri, 

Di = fii + (r-i^)flK, 
Ei — mi. 



(6) 



respectively, where K is the carrying capacity. The parameters 
where chosen to give modest positive growth rate (r = 1 .05 and 
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Figure 2: A: The average patch population size n at equilibrium, as a function 
of migration rate m from theory (solid line), and from simulations for K = 50 
(circles), A" = 30 (triangles), and K = \0 (diamonds). The parameters are N = 
1, 000, // = 1, and r = 1.05. B: The fraction of colonized patches as a function 
of migration rate m from theory (solid line), and from simulations (symbols). 
C: The average patch population size for an infinite number of habitats (dashed 
line) compared to simulations of finite number of habitats for K = 50 (circles, 
95% approximate confidence intervals). 



^ - 1) at low numbers. We will assume that births and deaths 
during migration can be ignored (v » 0), so that Eq. ^ applies. 
In order to investigate the accuracy of the theory, we will com- 
pare stochastic simulations of this dynamic to the solution of 
Eq. (|5]l, first with respect to the stationary states, for different m 
and K, and then to the full time evolution. 

4.1. Quasi- stationary states 

Fig. |2] shows the average patch population size, n, and the 
fraction of occupied patches as a function of the migration 
rate, from explicit simulations of the stochastic dynamic of 
N - 1, 000 patches and from numerical solution of Eq. (|5]l, 
for three different values of the carrying capacity K. The initial 
state for each simulation point was generated by distributing 
NK individuals randomly among the patches, so that the initial 
average number of individuals per patch is K. The populations 
were simulated to time 1 , 000, and the mean population size and 
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the fraction occupied patches were calculated as time-averages 
over the interval 200 < t < 1, 000. If the population went 
extinct before the end of the simulation interval, the mean pop- 
ulation size and fraction of occupied patches were taken to be 

zero. 

Hanski & Zhang ( 1993h analyzed a deterministic version 



of this model, where the local dynamic is replaced with 
the deterministic counterp art. Using the technique of 



Hanski & GvUenbera ( 119931) . they found that the dynamic con- 
verges to h = K and all patches occupied. When the migra- 
tion rate m is large, the average patch population size and the 
fraction of occupied patches are close to the values predicted 
by the deterministic local dynamic, but for small m the differ- 
ences are significant; especially, with stochastic local dynamic 
both patch population size and the fraction of occupied patches 
decrease with decreasing migration rate, and for each value of 
K there is a critical migration rate below which the population 
cannot sustain itself and becomes extinct. For all values of m, 
the stationary solution to Eq. (|5]l is in excellent agreement with 
the simulations. This illustrates the importance of accounting 
for stochasticity in the local dynamic, even though the global 
population changes deterministically. 

Fig. |2jZ! shows how the average patch population size typi- 
cally depends on the number of patches (see caption for param- 
eter values), compared to the theoretical value in a population 
with an infinite number of patches. The average is increasing 
with increasing number of patches, but also for small the rel- 
ative difference is not large. 




Figure 3: A: Comparison of the average patch population size in simulation 
trajectories of N = 10^ patches (black lines) to the solution of the nonUnear 
Master equation, Eq. (s), starting from the same initial distribution of patch 
population sizes as the simulations (grey lines), for three initial average patch 
population sizes (n = 5, n = 30, and n = 100). The parameters are = I, 
r = 1.05, K = 50, and m = 0.1. B: Time evolution of the fraction of empty 
patches, for the same trajectories as panel A. For h = 30 and n = 100 the 
fraction empty patches is approximately zero at all times. 



4.2. Transient dynamics 

An import advantage of the theory presented in this paper 
is that it not only gives the equilibrium points of the stochas- 
tic population dynamic, but the complete time evolution of the 
population. In this section, we compare the time evolution of 
the full stochastic simulations to the solution of Eq. (|5]l. 

Fig-EK shows the time evolution of the average patch popu- 
lation size from simulations of N - 1, 000 patches (black lines), 
with migration rate m - 0.1, starting from populations with 
three values of the average patch population size [n(0) - 5, 
n(0) = 30, and n(0) = 100]. The initial populations were 
produced by distributing Nh(Q) individuals randomly over the 
islands. Also shown are the corresponding solutions to the 
nonlinear Master equation, (Eq. |5] grey lines), starting from 
the same initial patch population size distribution as the cor- 
responding simulations. Simulations and theory are in close 
agreement, except some fluctuations in the simulations from the 
finite number of patches. These fluctuations become more pro- 
nounced around the stationary point of the dynamic. This can 
be expected since away from the equilibrium, relatively strong 
deterministic forces dominate the dynamic, but close to the 
equilibrium the dynamic is dominated by the stochastic com- 
ponents. 

Panel B of the same figure shows the time evolution of the 
fraction of empty patches for n - 5 (for the other values of 
n, the fraction of empty patches is approximately zero at all 
times), from simulations and the nonlinear Master equation. As 



before, there are small differences from the stochastic fluctua- 
tions in the simulation. The number of patches is initially in- 
creasing, before reaching a maximum and eventually decreas- 
ing towards zero. This is a marked difference from the deter- 
ministic theory, that predicts a steady decrease towards zero. 

Fig. ID shows the time evolution for a smaller migration rate, 
m = 0.01. This value is close to the critical migration rate for 
these parameters (c.f. Fig. As a consequence, the average 
patch population size at the stationary point is much lower than 
for m - Q.l 18 compared to ^ 42), and the fluctuations in the 
population are larger. Therefore, the simulation curves in Fig.|4] 
are averages over ten independent realizations of the stochastic 
dynamic, starting from the same initial distribution. 

As a final comparison, Fig.|5]shows the distribution of patch 
sizes in the simulations to the solution of the nonlinear Master 
equation at different times, from the simulations in Fig. [3] and 
Fig-El Because of the finite number of patches (A^ = 10^) there 
is significant spread around the prediction of Eq. (|5]l, but as can 
be seen in the right column, where the distribution is averaged 
over ten independent realizations, Eq. ^ predicts the expected 
distribution quite well. 

5. Discussion 

The model presented in this article makes several assump- 
tions that are idealizations of more realistic scenarios. Most of 
these are not necessary, but were introduced in order to sim- 
plify the description of the method. For example, it is assumed 
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Figure 4: A: Comparison of tlie average patch population size in simulation 
trajectories of N = 10^ patches (black lines) to the solution of the nonlinear 
Master equation, Eq. (5)> starting from the same initial distribution of patch 
population sizes as the simulations (grey lines), for three initial average patch 
population sizes (n = 5, n = 20, and h = 40). The parameters are fi = I, 
r = 1.05, K = 50, and m = 0.01. B: The time evolution of the fraction of empty 
patches. 
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Figure 5: Normalized histogram of patch population sizes from simulations 
(black lines) compared to the solution of the non-linear Master equation (grey 
lines). Left column: m = 0. 1, started from «(0) = 5, for t = \0, t = 50, and 
t = 100 (from the simulation in Fig.[3). Right column: m = 0.01 and «(0) = 40, 
for t = 2, t = 20, and t = 40 (from the simulation in Fig. |4). For increased 
legibility, the simulation histograms where binned in groups of two. 



that all habitats have the same size, availability of food etc., 
so that the dynamic is the same. In many cases, habitats are 
not equal but can vary significantly in terms of size and quality 
dHanski & GvUenbergl 1 19931: iHanskil Il994 . It is sti-aightfor- 
ward to extend the model presented in this paper to incorporate 
population structure, where the parameters governing the local 
patch dynamic may vary between patches. 

In addition, it is assumed that individuals migrate indepen- 
dently of others, although the rate of migration can have an arbi- 
trary dependence on the local patch population size. In biolog- 
ically more realistic situations, several individuals could leave 
and arrive together In addition, it may happen that both emi- 
gration and immigration leads to (or is caused by) a disruptio n 
of the local population, e.g. a conflict ( Johst & BrandlL[l997 ). 
It is possible to include such processes directly into the equa- 
tions of motion by modeling them with a probability transition 
matrix that connects the states before and after the event. 

Local colonization and extinction events are correctly cap- 
tured by Eq. (O, and if the dynamic is such that the popu- 
lation is driven to extinction in a deterministic manner, e.g. 
because deaths are systematically more frequent than births, 
these equations will show this as well. The risk of extinction 
of the whole population from demographic stochasticity cannot 
be calculated, however, since the derivations assume that the 
number of migrants in the dispersal pool is large. In order to 
study such events it is necessary to extend the calculation in the 
appendix to terms of the order ' . The resulting partial differ- 
ential equation is a diffusion approximation for the distribution 
of habitat population sizes, which can be used to approximately 



calculate the expected time to extinction of the whole popula- 
tion and more accurately estimate average population sizes. 

One of the main conceptual difficulties in modeling the pop- 
ulation dynamic in terms of population extinction and coloniza- 
tion events is in defining when an empty patch is 'co lonized' : 
Is it w hen the first individuals arrives at the path, as in (iKeeling , 



2002|), or should the patch be considered colonized first when 



the number of in dividuals is typical of a non-empty patch, as 



in jLande et al.V 11998,) ? The method in this paper avoids this 



issue completely, since it deals directly with the distribution of 
the number of individuals in the habitats. This is especially im- 
portant when habitats are small but many. In this case, there 
may not be a clear separation between colonized and empty 
patches, but rather a continuum of patch population sizes (c.f. 
Fig.EJ. 

The standard perception of the population dynamic of frag- 
mented populations has been that the local birth and death dy- 
namic is much more rapid than the rate of extinction and re- 
colonization of patches, so that the number of individuals in the 
occupied patches approaches a value determined by the number 
of occupied patches between each patch colonization or extinc- 
tion event. The derivations presented in this paper show that 
this view can be slightly misleading; rather than resulting from 
a separation of time scale between local dynamic and patch 
colonization and extinction, the effective population dynamic 
stems from a law of large numbers, where local fluctuations are 
buffered through a common migration pool so that the global 
dynamic remains smooth. Hence, the change in the fraction of 
occupied patches and in the average number of individuals on 
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occupied patches are coupled through their shared dependence 
on the rate of immigration. 

Finally, since the only critical assumption behind these re- 
sults is that the number of habitats is large, the theory is espe- 
cially suited to further elucidating the interplay of local stochas- 
tic dynamics and environmental factors in populations undergo- 
ing change. 
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A. Derivation of Eq. (|3) 

Our starting point is the time evolution of the probability p 
of having «, habitats with / individuals, for / = 0, 1, 2, . . . and M 
individuals in the dispersal pool, with the total count normal- 
ized as «, = A^. This gives an exact description of the full 
stochastic population dynamic. Taking into account the local 
dynamic in each patch, the migration process and the mortality 
due to migration, we obtain the Master equation 

d,p{no,n\,. . . , M) = v{M -H l)p(. . . ,M + 1) - vMp 

DO 

+ ^ \Tij{nj + l)p(. ..,ni-\,...,nj+\,...)- Tjitiipj 

i.j=0 



+ YjEi+i(ni+i + l)p(. - l,«,+i + 1,...,M- 1) 

1=0 

+ ^ a{M + \)-^ p(. . -H - 1, . . . ,M -H 1) 

!=1 

oo 

-J] +aM^ +vm)p (7) 



i=0 



for these probabilities. For brevity, we show only the vari- 
ables date differs from the left-hand side in these expressions. 
We now express the state in terms of the normalized frequen- 
cies fi - tii/N and the rate of immigration to a given patch, 
/ - aM/N. When the number of habitats is large, each event 
leads only to a small change in the probability p(/o,/i, ...,/) 
(each change is of the order 1/A^). Expanding p to order 
and ignoring terms of order 1 /N^ and higher, leads to the deter- 
ministic equations of motion for and /. 

5,p = J [r,, + 1 J |p - Ldjfi + i5,,p) - T,f,p 



N 



N 



N 



CO 

Collect terms in increasing order of we find that the coef- 
ficients of the leading term (proportional to A^) cancel, and the 
remaining terms can be written as 

oo r oo 

d,p = 2 ^'7 [dfMjp) - df,(fjP)] +d, (a + v)Ip -aj^ Ejrp 



ij=0 



+Yjdf, m + /) fip - ifi-ip - E^fi^ip] 



(9) 



where we have ignore all terms of order A^ ' or higher Hence, 
the probability obey a transport equation. Rearranging the 
terms, we obtain Eq. 
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